This tutorial is intended to show the basic outline of how OpenPNM works, and necessarily skips many of the more useful and powerful features of the package. So if you find yourself asking "why is this step so labor intensive" it's probably because this tutorial deliberately simplifies some features to provide a more smooth introduction. The second and third tutorials dive into the package more deeply, but those features are best appreciated once the basics are understood.
Learning Objectives
Python and Numpy Tutorials
Before diving into OpenPNM it is probably a good idea to become familar with Python and Numpy. The following resources should be helpful.
- OpenPNM is written in Python. One of the best guides to learning Python is the set of Tutorials available on the official Python website. The web is literally overrun with excellent Python tutorials owing to the popularity and importance of the language. The official Python website also provides a long list of resources
- For information on using Numpy, Scipy and generally doing scientific computing in Python checkout the Scipy lecture notes. The Scipy website also offers as solid introduction to using Numpy arrays.
- The Stackoverflow website is an incredible resource for all computing related questions, including simple usage of Python, Scipy and Numpy functions.
- For users more familiar with Matlab, there is a Matlab-Numpy cheat sheet that explains how to translate familiar Matlab commands to Numpy.
Before creating an OpenPNM simulation it is necessary to give a quick description of how data is stored in OpenPNM; after all, a significant part of OpenPNM is dedicated to data storage and handling.
OpenPNM employs 5 main objects which each store and manage a different type of information or data:
We will encounter each of these objects in action before the end of this tutorial.
Each of the above objects is a subclass of the Python dictionary or dict, which is a very general storage container that allows values to be accessed by a name using syntax like:
In [1]:
foo = dict() # Create an empty dict
foo['bar'] = 1 # Store an integer under the key 'bar'
print(foo['bar']) # Retrieve the integer stored in 'bar'
A detailed tutorial on dictionaries can be found here. The dict does not offer much functionality aside from basic storage of arbitrary objects, and it is meant to be extended. OpenPNM extends the dict to have functionality specifically suited for dealing with OpenPNM data.
All data are stored in arrays which can accessed using standard array syntax.
'throat.conns') that stores which pore indices are found on either end of a throat. This leads to an Nt-by-2 array. OpenPNM objects combine the above two levels of data storage, meaning they are dicts that are filled with Numpy arrays. OpenPNM enforces several rules to help maintain data consistency:
'pore.' or 'throat.'.'pore.diameter' or 'throat.length'.The following code snippets give examples of how all these pieces fit together using an empty network as an example:
In [2]:
import scipy as sp
import numpy as np
import openpnm as op
np.random.seed(10)
# Instantiate an empty network object with 10 pores and 10 throats
net = op.network.GenericNetwork(Np=10, Nt=10)
# Assign an Np-long array of ones
net['pore.foo'] = np.ones([net.Np, ])
# Assign an Np-long array of increasing ints
net['pore.bar'] = range(0, net.Np)
# The Python range iterator is converted to a proper Numpy array
print(type(net['pore.bar']))
In [3]:
net['pore.foo'][4] = 44.0 # Overwrite values in the array
print(net['pore.foo'][4]) # Retrieve values from the array
In [4]:
print(net['pore.foo'][2:6]) # Extract a slice of the array
In [5]:
print(net['pore.foo'][[2, 4, 6]]) # Extract specific locations
In [6]:
net['throat.foo'] = 2 # Assign a scalar
print(len(net['throat.foo'])) # The scalar values is converted to an Nt-long array
In [7]:
print(net['throat.foo'][4]) # The scalar value was placed into all locations
In [8]:
import openpnm as op
import scipy as sp
Next, generate a Network by choosing the Cubic class, then create an instance with the desired parameters:
In [9]:
pn = op.network.Cubic(shape=[4, 3, 1], spacing=0.0001)
The Network object stored in pn contains pores at the correct spatial positions and connections between the pores according the cubic topology.
shape argument specifies the number of pores in the [X, Y, Z] directions of the cube. Networks in OpenPNM are always 3D dimensional, meaning that a 2D or "flat" network is still 1 layer of pores "thick" so [X, Y, Z] = [20, 10, 1], thus pn in this tutorial is 2D which is easier for visualization.spacing argument controls the center-to-center distance between pores and it can be a scalar or vector (i.e. [0.0001, 0.0002, 0.0003]).The resulting network looks like: (This image was creating using Paraview, using the instructions given here)
In [10]:
print('The total number of pores on the network is:', pn.num_pores())
print('A short-cut to the total number of pores is:', pn.Np)
print('The total number of throats on the network is:', pn.num_throats())
print('A short-cut to the total number of throats is:', pn.Nt)
print('A list of all calculated properties is availble with:\n', pn.props())
In [11]:
print(pn.pores('left'))
The ability to retrieve pore indices is handy for querying pore properties, such as retrieving the pore coordinates of all pores on the 'left' face:
In [12]:
print(pn['pore.coords'][pn.pores('left')])
A list of all labels currently assigned to the network can be obtained with:
In [13]:
print(pn.labels())
In [14]:
geom = op.geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)
This statement contains three arguments:
network tells the Geometry object which Network it is associated with. There can be multiple networks defined in a given session, so all objects must be associated with a single network.pores and throats indicate the locations in the Network where this Geometry object will apply. In this tutorial geom applies to all pores and throats, but there are many cases where different regions of the network have different geometrical properties, so OpenPNM allows multiple Geometry objects to be created for managing the data in each region, but this will not be used in this tutorial.This freshly instantiated Geometry object (geom) contains no geometric properties as yet. For this tutorial we'll use the direct assignment of manually calculated values.
We'll start by assigning diameters to each pore from a random distribution, spanning 0 um to 100 um. The upper limit matches the spacing of the Network which was set to 0.0001 m (i.e. 100 um), so pore diameters exceeding 100 um might overlap with their neighbors. Using the Scipy rand function creates an array of random numbers between 0 and 0.0001 that is Np-long, meaning each pore is assigned a unique random number
In [15]:
geom['pore.diameter'] = np.random.rand(pn.Np)*0.0001 # Units of meters
We usually want the throat diameters to always be smaller than the two pores which it connects to maintain physical consistency. This requires understanding a little bit about how OpenPNM stores network topology. Consider the following:
In [16]:
P12 = pn['throat.conns'] # An Nt x 2 list of pores on the end of each throat
D12 = geom['pore.diameter'][P12] # An Nt x 2 list of pore diameters
Dt = np.amin(D12, axis=1) # An Nt x 1 list of the smaller pore from each pair
geom['throat.diameter'] = Dt
Let's dissect the above lines.
P12 is a direct copy of the Network's 'throat.conns' array, which contains the indices of the pore-pair connected by each throat.'pore.diameter' array, resulting in another Nt-by-2 array containing the diameters of the pores on each end of a throat.amin is used to find the minimum diameter of each pore-pair by specifying the axis argument as 1, and the resulting Nt-by-1 array is assigned to geom['throat.diameter'].'throat.conns' to index into a pore property array is commonly used in OpenPNM and you should have a second look at the above code to understand it fully.We must still specify the remaining geometrical properties of the pores and throats. Since we're creating a "Stick-and-Ball" geometry, the sizes are calculated from the geometrical equations for spheres and cylinders. For pore volumes, assume a sphere:
In [17]:
Rp = geom['pore.diameter']/2
geom['pore.volume'] = (4/3)*3.14159*(Rp)**3
The length of each throat is the center-to-center distance between pores, minus the radius of each of two neighboring pores.
In [18]:
C2C = 0.0001 # The center-to-center distance between pores
Rp12 = Rp[pn['throat.conns']]
geom['throat.length'] = C2C - np.sum(Rp12, axis=1)
The volume of each throat is found assuming a cylinder:
In [19]:
Rt = geom['throat.diameter']/2
Lt = geom['throat.length']
geom['throat.volume'] = 3.14159*(Rt)**2*Lt
The basic geometrical properties of the network are now defined. The Geometry class possesses a method called plot_histograms that produces a plot of the most pertinent geometrical properties. The following figure doesn't look very good since the network in this example has only 12 pores, but the utility of the plot for quick inspection is apparent.
In [20]:
water = op.phases.GenericPhase(network=pn)
Some notes on this line:
pn is passed as an argument because Phases must know to which Network they belong.pores and throats are NOT specified; this is because Phases are mobile and can exist anywhere or everywhere in the domain, so providing specific locations does not make sense. Algorithms for dynamically determining actual phase distributions are discussed later.Now it is necessary to fill this Phase object with the desired thermophysical properties. OpenPNM includes a framework for calculating thermophysical properties from models and correlations, but this is covered in :ref:intermediate_usage. For this tutorial, we'll use the basic approach of simply assigning static values as follows:
In [21]:
water['pore.temperature'] = 298.0
water['pore.viscosity'] = 0.001
We are still not ready to perform any simulations. The last step is to define the desired pore-scale physics models, which dictate how the phase and geometrical properties interact to give the transport parameters. A classic example of this is the Hagen-Poiseuille equation for fluid flow through a throat to predict the flow rate as a function of the pressure drop. The flow rate is proportional to the geometrical size of the throat (radius and length) as well as properties of the fluid (viscosity) and thus combines geometrical and thermophysical properties:
In [22]:
phys_water = op.physics.GenericPhysics(network=pn, phase=water, geometry=geom)
Network must be specifiedpores and throats where they apply, since this information is implied by the geometry argument which was already assigned to specific locations.We need to calculate the numerical values representing our chosen pore-scale physics. To continue with the Hagen-Poiseuille example lets calculate the hydraulic conductance of each throat in the network. The throat radius and length are easily accessed as:
In [23]:
R = geom['throat.diameter']/2
L = geom['throat.length']
The viscosity of the Phases was only defined in the pores; however, the hydraulic conductance must be calculated for each throat. There are several options, but to keep this tutorial simple we'll create a scalar value:
In [24]:
mu_w = 0.001
phys_water['throat.hydraulic_conductance'] = 3.14159*R**4/(8*mu_w*L)
Numpy arrays support vectorization, so since both L and R are arrays of Nt-length, their multiplication in this way results in another array that is also Nt-long.
Finally, it is now possible to run some useful simulations. The code below estimates the permeability through the network by applying a pressure gradient across and calculating the flux. This starts by creating a StokesFlow algorithm, which is pre-defined in OpenPNM:
In [25]:
alg = op.algorithms.StokesFlow(network=pn)
alg.setup(phase=water)
network argument.water, which dictates which pore-scale Physics properties to use (recall that phys_water was associated with water). This can be passed as an argument to the instantiation or to the setup function.Next the boundary conditions are applied using the set_boundary_conditions method on the Algorithm object. Let's apply a 1 atm pressure gradient between the left and right sides of the domain:
In [26]:
BC1_pores = pn.pores('front')
alg.set_value_BC(values=202650, pores=BC1_pores)
BC2_pores = pn.pores('back')
alg.set_value_BC(values=101325, pores=BC2_pores)
To actually run the algorithm use the run method:
In [27]:
alg.run()
This builds the coefficient matrix from the existing values of hydraulic conductance, and inverts the matrix to solve for pressure in each pore, and stores the results within the Algorithm's dictionary under 'pore.pressure'.
To determine the permeability coefficient, we must invoke Darcy's law: Q = KA/uL(Pin - Pout). Everything in this equation is known except for the volumetric flow rate Q. The StokesFlow algorithm possesses a rate method that calculates the rate of a quantity leaving a specified set of pores:
In [28]:
Q = alg.rate(pores=pn.pores('front'))
A = 0.0001*3*1 # Cross-sectional area for flow
L = 0.0001*4 # Length of flow path
del_P = 101325 # Specified pressure gradient
K = Q*mu_w*L/(A*del_P)
print(K)
The StokesFlow class was developed with permeability simulations in mind, so a specific method is available for determining the permeability coefficient that essentially applies the recipe from above. This method could struggle with non-uniform geometries though, so use with caution:
In [29]:
K = alg.calc_effective_permeability(domain_area=A, domain_length=L)
print(K)
The results ('pore.pressure') are held within the alg object and must be explicitly returned to the Phase object by the user if they wish to use these values in a subsequent calculation. The point of this data containment is to prevent unintentional overwriting of data. Each algorithm has a method called results which returns a dictionary of the pertinent simulation results, which can be added to the phase of interest using the update method.
In [30]:
water.update(alg.results())
Using Paraview for Visualization, the resulting pressure gradient across the network can be seen: